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ABSTRACT 

We implement the Elliptical Gauss-Laguerre (EGL) galaxy-shape measurement method pro- 
posed by Bernstein & Jarvis (2002) [BJ02] and quantify the shear recovery accuracy in weak lens- 
ing analysis. This method uses a deconvolution fitting scheme to remove the effects of the point- 
spread function (PSF). The test simulates > 10^ noisy galaxy images convolved with anisotropic 
PSFs, and attempts to recover an input shear. The tests are designed to be immune to shape 
noise, selection biases, and crowding. The systematic error in shear recovery is divided into two 
classes, calibration (multiplicative) and additive, with the latter arising from PSF anisotropy. 
At S/N > 50, the deconvolution method measures the galaxy shape and input shear to ~ 1% 
multiplicative accuracy, and suppresses > 99% of the PSF anisotropy. These systematic errors 
increase to ^ 4% for the worst conditions, with poorly resolved galaxies at S/N ~ 20. The 
EGL weak lensing analysis has the best demonstrated accuracy to date, sufficient for the next 
generation of weak lensing surveys. 

Subject headings: gravitational lensing — methods: data analysis — techniques: image processing 



1. Introduction 

Weak gravitational lensing, the shearing of 
galaxy images by gravitational bending of light, 
is an effective tool to probe the large-scale matter 
distribution of the universe. It is also a means to 
measure the cosniological parameters by compar- 
ing observation to numerical simulations of large 
scale structure growth (Bartelmann & Schneider 
2001). There are many weak lensing (WL) surveys 
underway to obtain the cosmological parameters 
to higher precision, and in particular to probe 
the evolution of the dark energy by observing 
its effects on the evolution of matter distribution 
(DLS^, CFHTLS2). 

The WL signal is very subtle, however; it is nec- 
essary to measure these small distortions (typical 
shear 7 ~ 1%) in the presence of optical distor- 
tions and the asymmetric point-spread-function 
(PSF) of real-life imaging. The level of system- 



atic error in the WL measurement methods are 
currently above the statistical accuracy expected 
from future wide and deep WL surveys (Pan- 
STARRS3, SNAP", LSST^, SKA^). Because there 
are no "standard shear" lenses on the sky, shear- 
measurement techniques are tested by applying 
them to artificial galaxy images and seeing if one 
can correctly extract a shear applied to the simu- 
lation. In most cases, the recovered shear can be 
written as 7out = '/niin + c. Departures from the 
ideal m = 1 we will term "calibration" or "mul- 
tiplicative" errors and quote as percentages. De- 
viations from the ideal c — can result from un- 
corrected asymmetries in the PSF and optics, and 
will be termed "additive errors" or "incomplete 
PSF suppression." 

Such tests of the most widely applied anal- 
ysis method (Kaiser, Squires, & Broadhurst 
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1995) [KSB], find m =0.8-0.9, but this coefficient 
is implementation dependent (Erben et al. 2001; 
Bacon et al. 2001), and depends upon the charac- 
teristics of the simulated galaxies. Hirata & Sel- 
jak (2003) [HS03] demonstrate that various PSF- 
correction methods can produce shear measure- 
ments miscalibrated by a few % to 20% or more. 
Heymans et al. (2005) [Shear TEsting Programme, 
(STEP)] present testing of many existing shear- 
measurement pipelines using a common ensemble 
of sheared simulated images. These methods show 
a median calibration error of 7%, although some 
(the BJ02 rounding kernel method, an implemen- 
tation of a KSB method, as well as the one de- 
scribed in this paper) show no calibration error, 
to within the <t„i w 1% noise level of the first 
STEP tests. Although the statistical accuracy in 
past surveys was comparable to the 7% systemat- 
ics, it is expected to be well below 1% in future 
surveys. Hence, understanding and eliminating 
the WL systematic errors require the most urgent 
attention today. 

In this paper, we implement the elliptical 
Gauss-Laguerre (EGL) deconvolution method as 
described in BJ02, and subject it to a series of 
tests designed to be more stringent than any pre- 
vious test of WL measurements. The deconvo- 
lution method is distinct from the Jarvis ct al. 
(2003) method, also described BJ02, in which 
the anisotropic PSF effects are removed using a 
"rounding kernel" instead. 

WL testing regimes are of two types: in end-to- 
end tests {e.g. STEP), one produces simulated sky 
images with a full population of stars and galax- 
ies, analyzes them with the same pipeline as one 
would real data, then checks the output shear for 
veracity. We perform here more of a dissection, in 
which we analyze the performance of the method 
one galaxy type at a time, and vary the param- 
eters of the galaxy and PSF images to determine 
which, if any, conditions cause the measurement 
to fail. While lacking the realism of an end-to-end 
test, this allows us to isolate and fix weaknesses. 
If we can demonstrate that the method succeeds 
under a set of conditions that will circumscribe 
those found on the real sky, then we can have con- 
fidence that our method is reliable, whereas end- 
to-end testing is reliable only to the extent that 
the simulated sky reproduces the characteristics 
of the real sky. 



Wc investigate here the performance of our 
EGL method across the range of noise levels, de- 
gree of resolution by the PSF, pixel sampling rates, 
galaxy ellipticity, and PSF ellipticity, using both 
highly symmetric and asymmetric galaxy shapes. 
We test not only the accuracy of shear recovery, 
but also the accuracy of the shear uncertainty es- 
timates. 

The EGL method is further elaborated in §2, 
while the implementation, GLFit, is detailed in 
§3. The shear accuracy test procedure is de- 
scribed in §4. The conditions under which the 
shape measurement succeeds, and the accuracy of 
its estimates of shear, are presented in §5. Pre- 
vious dissection tests include HS03 and Kuijken 
(2006). The former studies the performance of 
several methodologies on varied galaxy and PSF 
shapes/sizes in the absence of noise. The lat- 
ter study verified its "polar shapelet" method to 
better than 1% calibration accuracy. In §6 and 
§7 we conclude with comparisons to other shape- 
measurement methodologies and tests, and draw 
inferences for future surveys. 

2. Overview of Shape Measurement 

The task of this weak lensing methodology is 
to assign some shape to observed galaxy i, 
then to derive from the ensemble {e^} an estimate 
of the applied lensing shear 7. More precisely, 
a shape analysis can only determine the reduced 
shear 9 = 7/(1 — k), where k is the lens conver- 
gence. Following BJ02, we use distortion S to de- 
scribe the shear, where S = 2g/{l + g^) (S ~ 2g 
for g <C 1). In this paper, both the shear and the 
shapes are expressed as distortions; while in other 
WL literatures, shear is usually expressed as g. 

2.1. Defining Shapes using Shear 

Following BJ02, we will quantify the lensing by 
decomposing its magnification matrix M into a 
diagonal dilation matrix = e^l and a unit- 
determinant symmetric shear matrix Srj: 

M = D^Sr, = D^R^S^R_/3 (1) 

/ cos/3 -sin/3 \ , . 

^0 - [ sin/3 cos/3 ) 

= ( el^ ) 

V = (r?+,r/x) = (??cos2/3,r/sin2/3) (4) 
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where f3 is the direction of the shear axis, and 77 is 
a measure of shear. The "conformal shear" rj can 
be reparameterized as the distortion 6 = tanh r] or 
the reduced shear g = tanh 77/2. 

The shape must be assigned to an image of a 
galaxy with some surface-brightness distribution 
I{d). Initially we will ignore the effects of PSF 
convolution on the observed image. The B J02 def- 
inition of shape is to specify roundness criteria or 
"circularity tests," and Mx{I), that oper- 

ate on I to yield one scalar for each component 
of the shear — typically these are quadrupole mo- 
ments. The object is deemed circular (e = 0) if 
M+(7) = Mx (/) = 0. If the object is not circular, 
then we assign to the object the shape t] which 
yields the solutions 



M+[7(Sr,0)]=Mx[/(Sr,0)]=O, 



(5) 



i.e. we find the shear 77 that, when applied to 
the coordinate system 6, makes the image ap- 
pear circular in that coordinate system, and de- 
clare the galaxy shape to be this shear. Any cir- 
cularity test will do, as long as Equation (5) has a 
unique solution rj; in particular the matrix dM-fdrj 
must be non-singular. The shape e is defined by 
e = tanh 77, keeping the position angle f3. Defin- 
ing shape in this way with a suitable circularity 
test has the virtue that the effect of a lensing dis- 
tortion S upon the galaxy shape e is completely 
defined by the multiplication of shear matrices. In 
particular, the component-wise formulae for trans- 
formation of a shape under a shear must take the 
form given by Miralda-Escude (1991): 



e++S+ + {5x/S^)[l - Vl^]iex5+ - e+^x) . 



1 + e-S 



(6a) 



, _ ex + (5x + {S+/S^)[l - ^/^^](e+(5x - ex^+) 
1 + e-S 

(6b) 

We can take the limit of a weak shear 5+ <^ 1, 
<5x =0: 

eV = e+ + (1-6^)5+ (7a) 
e'x = ex-e+ex5+ (7b) 

BJ02 describe (§5) a scheme for optimally weight- 
ing and combining an ensemble of shapes to pro- 
duce an accurate estimate of the distortion S. 
This scheme is predicated on the assignment of 



shapes that transform under shear according to 
Equation (6). Hence to test the accuracy of our 
methodology in recovering weak lensing shear, we 
need only test that assigned shapes transform un- 
der shear as in Equations (7). Furthermore, the 
isotropy of the Universe guarantees that the Ci 
of an unlensed population will be uniformly dis- 
tributed in (3, hence we need only verify that Equa- 
tions (7) hold when averaged over an ensemble of 
galaxies with fixed unlensed |e| but random orien- 
tations: 

(e'+) = (e+) + (1 - = (1 - eV2)5(Sa) 

(e'x) = (e+) + (e+ex)<5+=0 (8b) 

Here the brackets refer to averaging over the pre- 
lensing orientation (3. The second order term in 
(5+ vanishes, so these equations arc valid to 0{S^). 
We refer to this as the ring test (Fig. 1), since we 
construct an ensemble of test galaxies which form 
a ring in the e plane, then shear them, measure 
their shapes, and take the mean. We also note that 
as a special case, we should obtain (e') = when 
there is no applied shear. When there is no PSF 
or the PSF is symmetric under 90° rotation, then 
this result holds for any measurement scheme that 
is symmetric under inversion or exchange of the x 
and y axes of images. But for an asymmetric PSF, 
this is a stringent test of the ability of the shape- 
measurement technique to remove the effects of 
the PSF from the galaxy shapes. 

2.2. Fitting to Basis Functions 

The circularity test described in BJ02 involves 
decomposing the pre-convolution surface bright- 
ness distribution of the galaxy, I{0), into the 
Gauss-Laguerre set of orthonormal basis functions 
in the plane. We consider here a general set of 
two dimensional functions {^pi{0)} that are com- 
plete (though not necessarily orthogonal) over the 
plane. Any set of complete functions can be trans- 
formed to a new complete set {ipf} via 



(0) ^ U^-'e) 



(9) 



Here E is any remapping of the sky plane that 
compounds a displacement xq, a shear r], and a 
dilation n. We refer to E as the basis ellipse of 
the function set: if the unit circle is an isophote of 
ipi, then {xajiJ-jT]} describe the center, size, and 
shape of an ellipse that traces the same isophote 
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for ^f. We will also use E to refer to this set 
of 5 parameters that defines an ellipse; it will be 
clear from context whether we are referring to the 
parameter vector or the coordinate transform. For 
any E we must have 



(10) 



5 = 





NX 


// 
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Fig. 1. — The ring test. In the ring test, objects 
that have ellipticity magnitude e = 



distributed evenly in /3 on a ring in the ellipticity 
plane (solid circle). In the absence of shear (5, the 
average of the measured shapes should be zero. In 
the presence of small shear 5+ , the distribution of 
objects in the ellipticity plane are displaced such 
that they are no longer arranged with their center 
on the origin. The average of all the shapes on the 
displaced ring should behave as Eq. (8), to the first 
order in (5+. 



for (at least) one vector h of coefficients. For the 
function set to be complete, the vectors must be 
infinite-dimensional, but a real measurement will 
model / with a finite subset of the basis functions. 
The model must be fit to the observed image plane 
which we denote with coordinates x. 

The action of the atmosphere, optics, and de- 
tector will operate on /(0) to produce an observed 
surface-brightness distribution loix) on the detec- 
tor plane. This observation operator C likely in- 
cludes convolution by the PSF and perhaps some 
distortion by the optics. We will assume that this 
operation is known and that it is linear over sur- 
face brightness, so that C{al\ -I- 6/2) = aC{l\) -I- 
hCil'i)-, where a and h are any two scalars. In this 
case, Equation (10) must also imply that 



(11) 

(12) 



The observed-plane brightness is sampled at the 
centers of the pixels Xp yielding measurements Ip 
with uncertainties tXp, assumed henceforth to be 
Gaussian. The model which maximizes the like- 
lihood of the decompositions (10) and (11), given 
the data is that which minimizes 



E 

pixels 



[lp-E^b.cl,f{Xp)Y 



(13) 



We call this a "forward fit" to the galaxy image, 
because we are positing a distribution b ■ tp^{6) 
of flux on the sky, then propagating this to the 
detector plane, where we compare the model to 
the observations. Our task will be to find the E 
for which the x^-minimizing vector b satisfies the 
circularity constraints.^ 



^Note this condition is subtly different from that imposed 
by Kuijken (2006): his ellipticity is that which minimizes 
if the b is always constrained to be circular. We do not 
yet understand the significance, if any, of this distinction. 
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For fixed E, the minimization is a linear 
least-squares problem over b with the usual an- 
alytic solution 

b = a-^-l3 (14) 

A = '£lpcl>f{xp)/al (15) 

p 

a,,- ^ ^cl,f{x,)cl>fix,)/al. (16) 

p 

If the ipi are orthogonal and the noise is station- 
ary and white and the sampling approaches the 
continuum limit and the PSF approaches a delta- 
function, then a will be diagonal (modulo some 
complex conjugation operations) and the solution 
is very simple. These conditions are not, how- 
ever, typically met by real data — in particular the 
sampling and PSF conditions — so it is not safe to 
use orthogonality to decompose the observed data 
(BJ02; Massey et al. 2004; Berry et al. 2004). 

For finite data, the solution must be done over 
some truncated basis set. We will assume that 
of the basis functions are being used, so that a, 
/3, and b have dimension N. The x^-minimizing 
b is uniquely defined, regardless of whether the ipi 
are complete or orthogonal, as long as a is not 
singular. 

To follow the BJ02 prescription we must de- 
fine circularity tests and iterate the shear in E 
until the two circularity components vanish. At 
the same time, we may wish to define null tests 
for the centroid xq and/or size fj, of the galaxy, so 
there will be some number K < 5 such statistics 
{K = 5 when there is a null test for each of the 
parameters of E). Any test which is linear over 
the true intensity I{6) may be expressed as a lin- 
ear operation on b. Hence the task of the fitter is 
to adjust K components of the E parameters until 
the x^-minimizing b satisfies 

M-b = 0. (17) 

M. is a, K X N matrix that defines our tests for 
matching the basis ellipse to the center, size, and 
shape of the pre-convolution galaxy. 

Once Equation (17) is satisfied for a chosen 
set of basis functions Vi and any choice of the 
circularity-test matrix M, we obtain a well-defined 
measurement of the shape e (or more generally the 
defining ellipse E) of the pre-convolution galaxy 
image. Since these shapes are defined by shear 



operations on the 6 plane, they should transform 
according to Equation (6) and be amenable to all 
the weighting and responsivity schemes of BJ02. 

Before describing how we solve for E, we sum- 
marize the conceptual and practical advantages 
of combining a transformation-based definition of 
shape with forward fitting of a pre-convolution 
model to the observed pixels. Note the fitting 
approach has been advocated by several authors 
(Fischer et al. 1997; BJ02; Massey et al. 2004; Kui- 
jken 2006). 

• The e values determined in this way trans- 
form in a well-known way under lensing 
shear. The response of a galaxy population 
to lensing shear can be calculated without 
recourse to empirical calibration factors be- 
yond the distribution of the e values them- 
selves. 

• The forward-fitting procedure can in princi- 
ple work with any sort of PSF, even those — 
such as the Airy function — that have diver- 
gent second moments. 

• The forward-fitting procedure properly han- 
dles pixelization and sampling by the detec- 
tor. Furthermore any aliasing ambiguities 
will properly propagate into the covariance 
matrix for b, and as described below can be 
propagated into uncertainties in E. 

• Pixels rendered useless by defects or cosmic 
rays are easily omitted from the measure- 
ment. 

• The method is easily extended to simultane- 
ous fitting of multiple images with distinct 
PSFs. The sum over pixels in Equation (14) 
is run over all pixels in all exposures of the 
galaxy. All exposures share the same b, ip, 
and E vectors, but have distinct observation 
operators C and hence different (jji. Infor- 
mation available in the best-seeing images is 
properly exploited. 

• The method is easily adapted to the analysis 
of uv-p\a,ne interferometric data rather than 
image-plane data. 

2.3. Iterating the Fit 

The desired ellipse parameters E = {xq, Uo, V+tVx} 
enter the least-squares fit non-linearly, hence some 
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iterative sehcmc (or Markov chain, cf. Bridle et 
al. (2002)) is required to determine the values that 
meet the constraint (17). For a chosen E, the de- 
termination of b has the rapid solution (14). If the 
current estimate Eq yields a 6o that does not meet 
the circularity condition, the Newton- Raphson it- 
eration would be 



En 



Eo + (5E, 6E 



dE 



(Mbo). 
(18) 

The derivative db/dE follows from noting the ef- 
fect of a small change SE to the basis of ip: 



E+5E 



(0) 



^f(5E-i0) 

OO 



OO 



^f+^^(x) 



k j=0 



Here is the generator for the transformation 
indicated by kth parameter of E — either transla- 
tion, dilation, or shear. These matrices are fixed 
by the choice of basis functions {ijji}. These al- 
terations to the basis-function values can be prop- 
agated through the solution (14) to give the per- 
turbation to b: 

6/3 = ^5EkGk0 

k 

Sac = Y,6Ek{Gk& + di^Gl) 

k 

b + Sb = (a-F<5a)-i(/3-F5/3) 
^ 6b = J2^^k[a~^Gk0-abo) 

k 

- a-^a'^Glbol .(19) 

The matrix db/dE is apparent from this last equa- 
tion. Here we have taken the generator matrices 
Gfc to each be A^x c»; the vector now must be, in 
general, augmented to infinite dimension, and we 
also take d to be oo x A^. Note that the parenthe- 
sized portion of the solution (19) would vanish if 
not for the distinction between the truncated and 
infinite-dimensional versions of a and /3, since the 
first N elements of /3 — abo are zero. Likewise 
we could set the initial a~^d^ of the final term 



to identity if not for the truncation, in which case 
the transformation would become the very simple 
Sb = (G^bo) • SE (with some abuse of notation 
here). Since we are using db/dE only to help us 
iterate the solution for E, we could use this simple 
approximation, or extend to some order beyond A'' 
using Equation (19). 

2.4. The Gauss-Laguerre Decomposition 

We use the Gauss-Laguerre decomposition in 

OUT shape measurements. These are the cigcnfunc- 
tions of the 2-dimensional quantum harmonic os- 
cillator, and are most compactly expressed as com- 
plex functions indexed by two integers p,q>0: 



2n 



(m = p—q > 0). 



(20) 



The elliptical-basis versions are taken to be 

^P^^{9) ^ e-2MV,p,(E-i0). (21) 

Note we have tweaked the normalizations so that 
the flux is 



f = J dHi{d) 



■ bpqSpq 
P>9>0 



(22) 



rather than making the functions orthonormal. 
The Gauss-Laguerre (GL) functions are still, how- 
ever, a complete and orthogonal set over the 
plane. The functions are equivalent to the "po- 
lar shapelets" of Massey et al. (2004). 

The advantages of the GL functions as our basis 
are: 

• The tpi are rapidly calculable via recursion 
relations. 

• Decompositions of typical ground-based 
PSFs and exponential-disk galaxies are rea- 
sonably compact. 

• There are obvious choices for the circularity 

tests: the centroiding condition is 6io ~ 0; 
the circularity condition is 620 = 0; and a 
size-matching condition is 611 = 0. The ma- 
trix M is very simple, i.e. it just restricts b 
to five relevant (real-valued) elements. 

• The generator matrices Gfc are sparse and 
simple. In particular, if the truncated b ex- 
tends to p-l-g' < N, then G(pg)(p'g/) vanishes 
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ioY p' + q' > N + 2. This means that the 
augmented a and /3 in Equation (19) can 
be calculated exactly simply by extending 
the sums in Equation (14) to order + 2 
(in p + q), just slightly beyond the order N 
needed to calculate b. 

• Matrices for finite transformations are cal- 
culable with fast recursion relations. 

• Convolution by a PSF is also expressible as 

a matrix operation over the b vector, if we 
also express the PSF as a GL decomposition 
6* in some basis ellipse E*. This means that 
the corresponding to a chosen PSF are 
rapidly calculable. 

BJ02 contains details on many of these proper- 
ties, and gives the generators (§6.3.2), the forms 
for the finite-transformation matrices (Appendix 
A), and the convolution matrices (Appendix B). 
We emphasize that other basis sets may be cho- 
sen which yield valid shape measurements, but we 
choose the GL decomposition for its efficiency and 
clarity. Similarly, one need not execute the oper- 
ation C{ippq) using the GL matrices; this is just a 
computational convenience in our case. 

2.5. Alternative Circularity Tests 

We follow BJ02 in using &20 = as our defini- 
tion of a circular object. It is possible that other 
GL coefficients carry shape information that could 
be used to produce lowcr-noisc shape estimates 
(Massey et al. 2004). In fact with the formalism 
used below (§2.6) to estimate shape uncertainties, 
one can derive for each object a linear circular- 
ity test M that is optimally sensitive to shear. We 
have implemented a scheme whereby such an opti- 
mally weighted combination of all the quadrupole 
coefficients {m = p — q = ±2) is derived for each 
galaxy. Our testing has not yet demonstrated 
any significant advantage to this scheme, how- 
ever, so we have not pursued it further for fear 
that such custom estimators might produce subtle 
noise-rectification biases. 

2.6. Error Estimates 

We may produce uncertainties on the ellipse pa- 
rameters Ek by propagating the pixel flux uncer- 
tainties. The covariance matrix for the b vector 



takes the usual form C** ^ a ^. The covariance 
matrix for the ellipse parameters then follows: 

H T 



fdb_ 
Vc/E 



fdb_ 



(23) 



We may also assign a detection significance to 
the object. In the simplest case, where there is no 
PSF to consider, the signal-to-noise ratio S/N for 
detection by the optimally-sized elliptical Gaus- 
sian filter is given by 



(24) 



Var(6oo) 

The maximization of this quantity with respect 
to the size, shape, and center of the elliptical- 
Gaussian filter generates conditions that are iden- 
tical to our default circularity tests. 

We define a significance in the more general 
case — such as when the PSF is non-trivial — to be 



Var(/)^ 



/ 



(25) 



The variance on the flux can be obtained by con- 
tracting C**. This estimate of the significance cor- 
responds to the signal-to-noisc on the galaxy de- 
tection that would be obtained by applying a filter 
to the deconvolved image that is matched to the 
shape, location, and radial profile of the galaxy. 
The maximum order of the sum in (25) determines 
the allowable complexity of the filter. 

3. Description of GLFit 

GLFit is the C++ implementation of the meth- 
ods described in the previous section. The code 
includes classes to represent the vectors, transfor- 
mations, and covariance matrices over quantities 
indexed by the integer pairs pq. To fit to a single 
galaxy, the code requires as input: 

• Image-format data for both the surface 
brightness Ip and the weight 1/(7^ (c/. 
Eq. (13)). The latter image should be set 
to zero at saturatc^d, defective, or contami- 
nated pixels. Data/weight image pairs from 
multiple exposures may be submitted for 
simultaneous fitting of a single object. It 
is assumed that each image has been flat- 
fielded so that a source of uniform surface 
brightness has equal value in all pixels. 
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• For each exposure, a map 6{x) from the 
pixel coordinate system into the world co- 
ordinate system. This is typically either the 

distortion map into the local tangent plane 
to the celestial sphere, or an identity map. 

• For each exposure, the local sky brightness 
and a photometric gain factor. 

• If we are executing a deconvolution fit, 
whereby we measure the shape of the galaxy 
before PSF convolution, then we require, for 
each exposure, a vector b* and basis ellipse 
E* describing the estimated PSF at the lo- 
cation of the galaxy. If we are executing a 
native fit, whereby we are trying to measure 
the observed shape of the galaxy, then no 
PSF information is required. 

• A choice for the initial order N {p+q < N) of 
the GL decomposition of the galaxy model. 

• Lastly an initial estimate of the size, loca- 
tion, and shape of the galaxy, e.g. the val- 

ited by detection software such as 
SExtractor (Bertin & Arnouts 1996). If we 
are executing a deconvolution fit, we need 
to know both the size (To = of the ob- 
ject as observed, and the size a^, of the PSF 
in which the object is observed. These are 
by default subtracted in quadrature to esti- 
mate the size ai of the intrinsic (pre-seeing) 
galaxy. 

The procedure for shape determination has the 
following steps: 

1. If we are executing a deconvolution fit, we 
choose the size of the "source basis" 
which is the basis ellipse for the ipf that will 
be used to model the pre- convolution image. 
We select 

e^^^ = a^^ = af + fpsFcrl (26) 

where af and are the sizes of the pre- 
seeing galaxy image and the PSF, respec- 
tively, and /psF is a prefactor chosen by de- 
fault to be 1.2. If there are multiple expo- 
sures, we take cr^^ to be the harmonic mean 
of all the exposures' PSF sizes. As discussed 
in BJ02 §6.2, we expect a choice of Gaussian 
basis a that is larger than both the intrinsic 



galaxy and the PSF to offer the most sensi- 
tive measure of pre-convolution quadrupole 
moments. During a deconvolution fit, /Xg is 
held fixed during all iterations and there is 
no constraint on size in the circularity test 
M. 

2. The initial guess for Eg is mapped from 
the world coordinate system into each ex- 
posure's pixel system. An elliptical mask is 
created for each exposure, defined by the no- 
contour of the elliptical Gaussian. Typically 
we select n « 4. Pixels outside the mask are 
not included in the fit. Larger masks reduce 
the possibility of degeneracy or non-positive 
models for high-order fits, but increase the 
chance of contamination by neighboring ob- 
jects. 

3. If any of the exposures do not fully contain 
the fitting mask, the Edge flag is set. If 
the initial centroid falls outside all the expo- 
sures, then the OutOf Bounds flag is set and 
the fit is terminated. 

4. The basis ellipse E* of the PSF is changed 
to have the same shape (not size) as the cur- 
rent estimate of the galaxy shape, by apply- 
ing a finite shear transformation matrix to 
b*. The convolution by the PSF is now ex- 
pressed as a matrix operation, so we have an 
observed-plane image model 

Io{x) = Y,4,f^[e{x)]Cijbj. (27) 

The observed-plane basis ellipse Eo has sec- 
ond moments that are the sum of the second 
moments of the chosen source-plane basis Eg 
and the PSF basis ellipse. The convolution 
matrix elements can be computed as in 
Appendix B of BJ02. We now have a model 
for the observed data which is linear over the 
expansion coefficients b. 

5. The linear solution for b is executed by stun- 
ming over all valid pixels within the masks 
on all exposures. If a singular value decom- 
position of the a matrix indicates that there 
are poorly constrained combinations of 6 el- 
ements, the order of the GL decomposition 
is reduced by one (to a minimum of A'^ = 2), 
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the ReducedOrder flag is set, and the solu- 
tion is attempted again. 

6. The matrix db/dE is calculated from the 
augmented d and /3 (which in fact were cal- 
culated in the previous step) as per Equa- 
tion (19). 

7. If the circularity test applied to b is not zero 
within a specified tolerance, then an itera- 
tion increment SE is calculated according to 
Equation (18). If a singular- value decompo- 
sition of db/dE indicates a nearly-degenerate 
matrix, we know that there is no sensible es- 
timate for the next step. In other words, at 
least one of our circularity tests is currently 
ill- defined, for example no small change of 
basis makes the object appear more circular. 
This should, for example, be the case for a 
point source, but also arises for some very 
irregular or noisy galaxies. The behavior of 
the next step depends upon the fit type: In 
a native fit, we set the FrozeDilation flag, 
set fi back to its initial value, and restart the 
fitting process, this time holding /x fixed. If 
we are executing a deconvolution fit, we in- 
crement /psF by 0.5, set the RaisedMu flag, 
and restart the fitting process, since fits to 
a larger basis tend to be more stable (but 
less sensitive). If FrozeDilation has been 
set (for a native fit) or /psF > 3 already 
(for a deconvolution fit), and db/dE is still 
near-degenerate, then our last resort is to set 
the FrozeShear flag, reset r] to its starting 
value, and restart the fit process while hold- 
ing the shape fixed. 

8. The next iteration of the fitting process con- 
tinues at step 4. If the circularity test is 
not satisfied within a chosen number of iter- 
ations, the DidNotConverge flag is set and 
the fit is terminated. If any other ma- 
trix singularities are encountered, or if there 
are too few pixels to conduct the fit, the 
Singularity flag is set, and the fit termi- 
nates. 

9. Once the circularity tests are satisfied within 
a desired tolerance, the resultant E and b for 
the galaxy are available for the user. The a 
matrix provides the covariance matrix for b. 



The errors on the ccntroid, size, and elliptic- 
ity (or whichever of these were free to vary) 
are calculated using Equation (23), and a 
detection significance can be reported using 
Equation (25). 

In most cases the solution for E is executed twice: 
first, a coarse fit at order N = 2 is attempted. 
If this fails, the CoarseFailure flag is set. If it 
succeeds, the coarse solution is used as the starting 
point for a fit to the full requested N. 

Another subtlety is that the pixel mask is not 
changed between iterations — the mask maintains 
the shape and size specified by the initial guess (or 
the result of the coarse fit). This should not bias 
the fit toward the initial guess, however, because 
we are fitting to the data within the mask rather 
than executing a sum of moments over pixels 
within the mask. If the size /x of the object changes 
dramatically during the coarse fit [e.g. due to a 
very poor initial estimate from SExtractor), then 
we do resize the masks and start over. If at any 
time the mask on an image grows to encompass too 
large a number of pixels, we set the TooLarge flag 
and quit. This prevents the program from getting 
hung up calculating pixel sums for nonsensically 
large objects. 

The algorithm is considered to have con- 
verged to a valid shape measurement if it com- 
pletes without setting of the flags FrozeShear, 
DidNotConverge, Singularity, OutOf Bounds, 
TooLarge , or CoarseFailure. 

4. Test Procedure 

In our tests, postage-stamp images arc created 
of individual galaxies, as convolved with a known 
PSF, pixelized, and given noise. We then extract 
shape and shear measurements from ensembles of 
such postage-stamp tests. In the real sky and in 
end-to-end tests which draw (pre-lensing) galaxy 
shapes at random, the uncertainties in the output 
shear are usually dominated by the finite variance 
of the mean pre-shearing shape (shape noise). We 
construct our ensembles with zero mean shape to 
eliminate shape noise, leaving our testing precision 
dominated by the shot noise in the images. 

Because we provide the shape-measurement al- 
gorithms with an approximate location for every 
simulated galaxy, our tests do not include the 
eff'ects of possible selection biases (Kaiser 2000; 
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BJ02; HS03). Similarly the PSF is always pre- 
sumed to be known exactly, so we do not test for 
errors that would result from PSF measurement 
or interpolation errors (van Wacrbckc ct al. 2005; 
Jain et al. 2006). Effects of galaxy crowding or 
overlap, detector nonlinearities, intrinsic galaxy- 
shape correlations, or redshift determination er- 
rors (Ma et al. 2006) are not examined here. In- 
stead we concentrate on the sources of system- 
atic error that arise specifically from (1) the shape 
measurement, (2) removal of the PSF effects, and 
(3) the shear estimation, in images with finite 
noise and sampling. 

4.1. Test Conditions 

In order to test our WL analysis methodology, 
we need to find the conditions imder which GLFit 
works, and then test for any biasing in the shape 
measurement or shear estimation method. The 
following is the list of the test conditions or choices 
made: 

• Native fit or deconvolution fit. All measure- 
ments are tested for cases without (native 
fit) and with (deconvolution fit) the presence 
of a PSF. In a WL analysis of real sky data, 
native fits are used in (1) characterizing the 
PSF from stellar images, and (2) measuring 
the shape of the PSF-smcarcd galaxy, whose 
value is used to produce an initial estimate 
of the intrinsic galaxy shape. 

• Galaxy size with respect to pixel or PSF. 

As the pixel sampling rate decreases, it be- 
comes harder to obtain shape information 
about the galaxy from the pixelized image, 
or for the GLFit to converge. For the native 
fits, the galaxy size is defined with respect to 
the pixel; for the deconvolution fits, to the 
PSF. 

• Galaxy detection significance v, or equiva- 
lently S/N. If the background noise dom- 
inates the signal from the galaxy, GLFit is 
also unlikely to converge. 

• Galaxy shape e. The galaxies are sheared to 

various shapes, ranging from |e| = to 0.9 
(= axis ratio 4.35 : 1). In some cases, the 
major axis orientation (parallel or diagonal 
to the pixel axis) is also examined. 



• Galaxy type. We choose a "symmetric" or 
"asymmetric" galaxy model, with the latter 
being a more stringent test of the method- 
ology. The galaxy models arc described fur- 
ther in §4.2, and the shear accuracy is exam- 
ined by the ring test, described in §4.3. 

For the deconvolution fit, we have the following 
additional choices for the PSF: 

• PSF type. We choose an Airy disk PSF as 

the worst possible case for a PSF. An Airy 
disk offers particular difficulties for shape- 
measurement algorithms, such as its non- 
trivial morphology and divergent second mo- 
ments. 

• PSF shape. The Airy disk is either circular 
or elliptical (sheared at epsF = 0.1). We 
look for any leakage of the PSF shape into 
the galaxy shape measurement or the shear 
estimate. 

In this paper, the PSF is always known; we do 
not test for errors due to incomplete knowledge of 
the PSF. 

4.2. Pixelized Galaxy Models 

Two types of galaxies, symmetric and asym- 
metric, are used to test the shape measurement 
and shear recovery. The galaxy models, along with 
any shear, rotation, and dilation, are described as 
an analytic function over the pixel coordinates. 

The "symmetric" galaxies have either a Gaus- 
sian or exponential circularly symmetric radial 
profile that is then sheared. All isophotes are sim- 
ilar ellipses, so these galaxies have unambiguous 
shapes, which allows direct comparison between 
the input and measured shape. The elliptical sym- 
metry of these galaxies could be masking errors in 
the shape-measurement methodology by causing 
fortuitous cancellations of errors. In other words, 
all of their EGL expansion coefficients with m 
are exactly zero, so any biases that couple to these 
coefficients will not be tested by the symmetric 
galaxies. 

The "asymmetric" galaxies are designed so 
that they have no inversion symmetry, and their 
isophotal shape varies with radius; i.e., there is 
no unambiguously defined shape. The ring test is 
therefore needed to see if input shear is accurately 
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recovered. These are constructed as a combi- 
nation of exponential and deVaucouleur ellipses, 
each with different flux, size, and shape, as well 
slight offset between the two ccntroids. The 
asymmetric galaxies offer a more stringent test 
of shear measurement, since we have broken any 
symmetries that could mask errors in the shape 
measurement. 

A pixelized postage stamp image is created 
from the analytic functions that describe the 
galaxy and PSF models at the intended sampling 
rate. The centroid of the analytic functions is ran- 
domized within the pixel over multiple image re- 
alizations (pixel phase randomization). The PSF 
convolution is done either analytically if possi- 
ble, or via FFT of noiseless pixelized images. The 
Poisson photon noise is then added to produce the 
final image to produce the desired S/N for the ob- 
ject. We then nm GLFit for shape measurement 
with starting parameters randomly displaced from 
the true values. 

4.3. The Ring Test 

The ring test evaluates the mean shape — the 
weighted average of the galaxy shape components — 
with a simplified galaxy shape distribution, as seen 
in Figure 1. When a shear S is uniformly applied 
to the ensemble, the mean shape (e) should equal 
nS, where 7^ = 1 - ^ (^q- (8)) is called the rc- 
sponsivity. This relation then gives an estimate of 
the shear from the mean shape. 

The value of the ring test is two-fold: the ex- 
pected signal is known, and the shape noise is 
absent. The shape noise in WL is the statisti- 
cal noise due to the random distribution of shapes 
with finite magnitudes; the uniform distribution of 
shapes in the ring test eliminates this noise. The 
lack of shape noise and the known expected value 
then allow systematic errors to be quantified. 

5. Results 

5.1. GLFit Convergence 

We run GLFit on objects with different size and 

detection significance to understand the condi- 
tions under which GLFit converges reliably, which 
we define to be > 99% rate of successful shape 
measurement. Shear estimate accuracy tests are 
subsequently performed only in > 99% conver- 



gence conditions, to preclude any selection biases 
that could be generated by the non-convergence of 
GLFit. 

5.1.1. Native Fit Convergence 

Figure 2 shows the 99% convergence contour for 
native fits on exponential-profile symmetric galax- 
ies of various ellipticities. The upper right hand 
region delineated by the contour indicates where 
the fit has > 99% convergence. The horizontal 
axis indicates the sampling rate, and is expressed 
in terms of the "minor axis width" (MAW), de- 
fined below. The vertical axis incidates the detec- 
tion significance v of the object. 

We define the minor axis width for a Gaussian 
object as its full width at half maximum (FWHM) 
along the minor axis direction. For non-Gaussian 
objects, the "Gaussian size" a is defined from the 
size-matching condition fen = (§2.4); an expo- 
nential profile /(r) = /oCxp(— r/ro) has a Gaus- 
sian size cr ~ 1.16ro. The Gaussian size along the 
minor axis is (Jb = crexp(— 77/2), where r] is the 
conformal shear (§2.1), so the minor axis width is 
2.35(76. 

Figure 2(a) clearly shows that for a native fit 
to converge reliably, S/N must be > 10 and its 
minor axis width be > 2.8 pixels. The convergence 
contour for Gaussian radial profiles is found to be 
identical to the exponential-profile contour. The 
fact that the minor axis dimension is the limiting 
factor for convergence implies a shear-dependent 
selection, since <Jb does not remain constant under 
a shear operation. 

The required sampling rate depends, however, 
on the orientation of the minor axis. Figure 2(a) 
shows the convergence contour when the ellipse is 
elongated along the pixel axis, where Figure 2(b) 
shows the case when the elongation is along the 
pixel diagonal. In the second case, the required 
sampling rate appears to decrease with increasing 
e, until it settles at 2.0 pixel for e > 0.6. Figure 3 
illustrates that for objects elongated along the 
pixel diagonal, the effective sampling rate along 
the minor axis becomes l/\/2 of the pixel spacing. 

The GL decomposition of the PSF is deter- 
mined via a native fit. Hence, the images must 
sample the PSF minor axis width by more than 
2.8 pixels for the GLFit to be applied successfully 
to individual stars. Note that it is possible to in- 
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crease the sampling rate by dithering for a CCD 
of any pixel size. 



5.1.2. Deconvolution Fit Convergence 

For the deconvolution test, we convolve the 
model galaxies with a circular Airy-function PSF, 
whose characteristic width is ctpsf — 0.21 X/D. 
The model galaxies in this test are symmetric el- 
liptical objects of exponential radial profile with 
minor axis Gaussian width tTf, before PSF convolu- 
tion. The pixel size is irrelevant to the convergence 
rate, as long as crpsF > 2.8 pixels (c/. §5.1.1). We 
find that the convergence rate of GLFit is inde- 
pendent of the orientation of the major axis under 
these conditions. 

Figure 4 shows the 99% convergence contour 
for deconvolution fits. The horizontal axis is the 
minor-axis resolution 



rb = CTb/fpsF- 



(28) 



We find that high-e galaxies are less likely to con- 
verge at low resolution than more circular galaxies. 

In analyses of the real sky, the dimmest and 
smallest galaxies are most numerous, so we draw 
our attention to the lower left corner of the con- 
tour in Figure 4. GLFit will converge reliably for 
objects with rf, > 1.8 down to > 20. To obtain a 
reliable shape measurement of marginally resolved 
objects (rf, > 1.2) at all shapes e, the significance 
must be > 40. 

5.2. Accuracy for Symmetric Galaxies 

Having determined where GLFit is consistently 
successful in measuring a shape, we now examine 
the accuracy of the shape measurement and error 
estimate that result. 

When an elliptical object is used as an input, 
the input shape is well defined, and hence so is the 
measurement error. In this section we use sym- 
metric elliptical galaxies with exponential radial 
profile that have an unambiguous shape vyinput- 
The shape measurement error is expressed as a 
conformal shear A77 = (— Tyinput) © Vmeas, where ® 
is the shear addition operator, equivalent to shear 
matrix multiplication: 



Sr,^R=S^^S^^, (29) 



where R is the unique rotation matrix that al- 
lows Srj^ to be symmetric (BJ02 §2.2). The quan- 
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Fig. 2. — Contours of 99% successful GLFit 
native-fit convergence over pixel sampling (hori- 
zontal axis) and object S/N (vertical axis). The 
pixel sampling is defined as the FWHM along the 
minor axis direction, or the "minor axis width" . 
The shaded region is where the native-fit conver- 
gence rate is < 99%. (a) Convergence when major 
axis is along the pixel axis, (b) Convergence when 
major axis is along the pixel diagonal. The shift 
in the contour in (b) is due to an effectively new 
pixel spacing (c/. Fig. 3). The circles in each plot 
indicate the parameters at which the shape mea- 
surement accuracy tests have been performed (c/. 
Fig. 5). 
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Fig. 3. — Relation of the minimum mcasureable 
object size (Fig. 2) to the pixel lattice for native 
GLFit. The circle and the ellipse have the mini- 
mum minor-axis width for reliable shape measure- 
ment in native fits. For the ellipse on the right 
(cx = 0.6), the sampling spacing along the mi- 
nor axis, a', is effectively l/\/2 of the pixel lattice 
spacing a. The sampling spacing is still a for el- 
lipses stretched along the pixel axis. We note that 
the the object center is placed at random locations 
within the pixel in our simulations. 



tity Ary is the amount of shear required to bring 
the input shape to the measured shape, hence de- 
scribes a shear bias regardless of the intrinsic ob- 
ject shape. Each data point is an average over 
N = 100, 000 measurements of independent re- 
alizations of pixel phase and photon noise. We 
then plot the averaged error components (A77+) or 
(A?7x ) as the measurement accuracy. The uncer- 
tainty in the mean + component is then cr^,,^^^ = 

GLFit produces error estimates (t^_^ , ct^^ for ev- 
ery shape measurement from the covariance ma- 
trix (§2.6). We test for the accuracy of 
this estimation by examining the averaged ratios 
((A?7+)^/ct^^) and ((Aryx )^/(Tj^^ ), where the aver- 
age is over shape measurement trials. 

5.2.1. Native Fit Errors 

Each panel in Figure 5 show the average error 
(Aryi) (i = x) as a function of the input shape 
e of the symmetric exponential galaxies. We first 
examine the accuracy of native fits to galaxies with 
no PSF convolution. The panels correspond the 
sampling rates and detection significances marked 
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Fig. 4. — Contour for 99% deconvolution-fit shape 
measurement convergence over resolution (hori- 
zontal axis) and object S/N (vertical axis), for an 
Airy PSF convolution. The shaded region is where 
the deconvolution-fit convergence rate is < 99%. 
The resolution is expressed as = (Jf,/crpsF, 
where upsF and Ob are the characteristic Gaus- 
sian widths for the PSF and along the galaxy mi- 
nor axis, respectively. The contour is independent 
of the orientation of the major axis. The circles 
indicate the parameters at which the shape mea- 
surement accuracy tests have been performed (c/. 
Fig. 8). 

by the points inside the 99% convergence contour 
in Figure 2. 

Those that are well into the > 99% conver- 
gence region (the inner three panels) show no sys- 
tematic bias in the shape measurement within the 
uncertainty tT(^.)< 0.0005. For those on the 99% 
convergence contour (left-most and bottom pan- 
els), shape measurements show some systematic 
errors, but the worst systematic is still within 
A?7i/ei ~ 0.4%. The shear error A7/7 (— m — 1 if 
c = 0) that would result from such measurement 
error is 

A7 _ A77, 
— / 7 

7 

assuming /S.5/5 ~ A7/7, where 



(30) 



/ 



1 



(e2)/2 



1 



(e2)/2 



(l-(e'» 



(31) 



for a constant l^rji/ci. The function / is always 
< 1, and rapidly decreases to zero for 0.6 <e < 
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Fig. 5. — The error in native fit shape measure- 
ment as a function of the input shape e. The 
"shape error" A77 is the amount of shear necessary 
to bring the measured shape to the input shape; 
Arji is the ith component of the error. The plots in 
the same column have the same minor axis width 
(MAW= 2.350-;,), whereas the plots in the same 
row have the same significance v] each panel cor- 
respond to the set of {v, MAW) parameters indi- 
cated by the circles in Figs. 2(a) and (b). The 
slanted, dotted line in two of the panels indicate 
shape error of Aiji/a — —0.4%, where i = +,x. 

1. Hence a shape measurement systematic of 
Ar]i/ei ~ 0.4% approximately corresponds to a 
shear calibration error of A7/7 ~ 0.4%. 

Figure 6 shows the root-mean-square (RMS) of 
the actual to estimated error ratios, ^ {Arjf / ct^ . ) , 
where ct,,. is the error estimate in the ith com- 
ponent as calculated by GLFit. Overall, we see 
that the error estimate does fairly well. There is 
a slight tendency to overestimate the error as 
becomes low; the ct^^ bias has no dependence on 
the pixel sampling. The v dependence is approxi- 
mately ^(aJ/^c. 1-M. 

5.2.2. Deconvolution Fit Errors 

Figure 7 shows the accuracy of a deconvolution 
fit when a symmetric, exponential galaxy is con- 
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Fig. 6. — The RMS average of the ratio of mea- 
sured error to estimated error. The error estimate 
arii is obtained from the fitting process (c/. §2.6). 
Each panel correspond to the set of parameters 
indicated by the circles in Figs. 2(a) and (b). 



volved with an Airy PSF. We first discuss a high- 
S/N case, v = 100. The measurements were done 
at various minor axis resolutions, where = 0.5, 
1, 2, and 5 each correspond to the columns from 
left to right. 

The panels in the top two rows show the shape 
error (Ar/i) as a function of the input galaxy el- 
lipticity. In the top row, the Airy PSF is cir- 
cular (epsF = 0.0), while in the middle row the 
Airy PSF is anisotropic, with an ellipticity of 
epsF = 0.1. The PSF and galaxy ellipticity both 
are elongated along e+ ; we find that the results do 
not change even if we add a non-zero ex < 0.15 
component to either the PSF or the galaxy. The 
major axis orientation with respect to the pixel 
axis is irrelevant for deconvolution fits (c/. §5.1.2). 

The bottom row plots the difference in average 
measured shape between epsF = 0.1 and epsp = 
0.0. If the PSF anisotropy is fully suppressed in 
the deconvolved galaxy-shape measurement, this 
difference should be zero. 

In general. Figure 7 shows the following two 
systematic shape measurement errors: (1) Regard- 
less of the PSF anisotropy, the measurement er- 
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Fig. 7. — The error in deconvolution shape mea- 
surement at v = 100, at resolutions = 0.5, 1, 
2, and 5. The plots in the same column have the 
same resolution. The top row is convolved with 
a circular Airy PSF; middle row with an ellipti- 
cal Airy PSF (epsp = 0.1); the bottom row is the 
difference between the two. In the top two rows, 
the slanted, dotted lines indicate shape error of 
Ar7+/e+ = ±1%, where in the bottom row, the 
dotted lines indicate = ±0.001, or 1% of the 

PSF anisotropy. For rf, = 5, the systematic specif- 
ically due to the GL decomposition of the Airy 
PSF is removed by intentionally blurring both the 
PSF and the observed image with a Gaussian that 
is a fraction of the size of the galaxy. 








Fig. 8. — The error in deconvolution shape mea- 
surement for 1/ < 50. The Airy PSFs are isotropic 
(epsF = 0). The plots in the same column have the 
same resolution rf,, whereas the plots in the same 
row have the same significance each panel cor- 
respond to the set of (i^, rt,) parameters indicated 
by the circles in Fig. 4. These plots are similar to 
those on the top row of Fig. 7. The slanted, dotted 
lines indicate shape error of Arj^/e^ = ±1%. 
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Fig. 9. — The PSF anisotropy bias in deconvolu- 
tion shape measurement for < 50. Each panel 
correspond to the set of (v, rt) parameters indi- 
cated by the circles in Fig. 4. These plots are 
similar to those on the bottom row of Fig. 7. The 
dotted lines indicate (A?7+) ~ ±0.001, or 1% of 
the PSF anisotropy. 
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ror exhibit a slight systematic when the resolu- 
tion is low (rf, < 2). There is a constant positive 
slope in (A77+) between < e < 0.6. (2) For 
e > 0.6, the slope turns around, eventually un- 
derestimating the shear at e = 0.9. In the region 
described by (1), the shear measurement errors 
are approximately (Ayyi) < O.Ole^, corresponding 
to A'y/j< 1%. The well-resolved galaxies {ri, > 5) 
exhibit biases in shape measurement, which are a 
consequence of a peculiarity of the Airy function 
PSF. We will explain this below, and describe a 
remedy which yields the improved performance of 
the points marked "blur." 

Despite these measurement systematics, how- 
ever, the difference between the epsp = 0.1 and 
0.0 measurements are less than a shear of 0.001, 
equivalent to < 1% of the PSF anisotropy, except 
for the poorly resolved {vb = 0.5) and highly el- 
liptical galaxies. Hence our GLFit implementation 
has a PSF anisotropy bias suppression of > 99% 
at 1/ = 100, except for the poorly-resolved case. 

We now examine deconvolution fits with galax- 
ies of lower S/N, using the combinations of reso- 
lution and significance marked by the circles sam- 
pling the > 99% convergence region in Figure 4. 
Figure 8 shows the measurement error (A?7_|_) for 
a circular PSF, and Figure 9 shows the shape 
bias induced by a PSF ellipticity of epsF = 0.1. 
The content of these figures are similar to those 
on the top and bottom row of Figure 7, respec- 
tively. Evaluated at input e = 0, the initial slope 
in (A77_|_) and the value of the PSF bias increase 
as '--^ l/i^ as J' decreases. In both cases, there is a 
turn-around in the magnitude of the systematic, 
with a zero crossing in the region 0.6 ^e. A rough 
estimate for the shear error A7/7 due to these 
shape errors (assuming a flat e distribution) is 1- 
2%. 

Figure 10 plots the RMS error ratio of the ac- 
tual to the estimated ^ {Arj"^ / a'^^) for deconvo- 
lution fits, and is analogous to Figure 6. For the 
deconvolution, the error estimate ct^^ does fairly 

well for > 50, while ^{Arj^/a^^) ~ 1 — | when 
I' < 50. 



Fig. 10.— The RMS ratio of actual {A7J+) to es- 
timated (o'ri^) errors for deconvolution fits. Each 
panel correspond to the set of {i/, r^) parameters 
indicated by the circles in Fig. 4. 



5.2.3. The Airy Fix 

The open-circle points in the rightmost column 
of panels in Figure 7 mark a tendency for the de- 
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convolution procedure to underestimate the ellip- 
ticity of large objects. This can be traced to the 
fact that the Airy function is rather poorly de- 
scribed by a Gauss-Laguerre expansion. In our 
default procedure, the size of the basis set of the 
EGL expansion of the PSF is made similar to the 
size of the PSF itself; in this case the Airy function 
is poorly modelled at radii ^ X/D. For example 
the Airy function has divergent second radial mo- 
ment, while any EGL expansion at finite order has 
a finite second radial moment, indicating that the 
EGL expansion has failed to capture the large-r 
behavior of the Airy function. It takes GL order 
A^GL = p + q = 8 and 16 to describe the first and 
second Airy ring, respectively; our simulations use 
GL order 12 to describe the PSF. 

An equivalent statement is that, in the Fourier 
domain, an EGL expansion fails to properly de- 
scribe the small-fc behavior of the Airy function. 
This is not surprising since the Airy function has 
a cusp at fc = (which is equivalent to having in- 
finite second moment in real space). The cusp in 
turn results from the sharp edge of the illumina- 
tion function of an unapodized circular mirror. 

The poor description of the Airy function at 
large r or small k becomes important when trying 
to deconvolve images of well-resolved galaxies, be- 
cause the shape information for large galaxies is 
carried at large r and small k relative to the PSF. 
The finite EGL expansion of the Airy PSF under- 
estimates its circularizing effect on large galaxies, 
hence the deconvolved shapes are too round for 
large galaxies. Marginally-resolved galaxies don't 
have this difficulty because they carry their ellip- 
ticity information in the part of fc-space where the 
EGL expansion is a good match to the Airy func- 
tion. 

We sec that this "Airy failure" is not an intrin- 
sic difficulty of the shape-measurement methodol- 
ogy, but rather stems from a poor model of the 
PSF. This leads us to consider several possible so- 
lutions: 

• Describe the PSF using a more appropriate 
function set than the Gauss-Laguerre expan- 
sion, if one expects nearly diffraction-limited 
images. This would complicate numerical 
calculation of the operator C from §2.2, but 
should fully restore the numerical accuracy 
for any PSF. 



• When fitting well-resolved galaxies, choose 
a basis for the EGL expansion of the PSF 
which matches the size of the galaxy rather 
than the size of the PSF. This puts the fit- 
ting freedom of the EGL expansion in the 
range of fc-space that is more important for 
the problem at hand; the core of the Airy 
function will be poorly fit but the wings will 
be better. 

• Change the PSF. 

In these tests we implement the last of these three 
options, by blurring the postage-stamp and PSF 
images with a Gaussian that is a 1/4 of the size 
of the galaxy, and then apply GLFit deconvolu- 
tion using the smeared image pair. The blurring 
decreases the size mismatch between the PSF and 
galaxy GL bases, which improves the accuracy of 
the GL expansion at the size scale relevant to the 
large galaxy. The blurring very slightly increases 
the scatter of the shape measurement, but greatly 
reduces the bias, as seen in the figure. 

To summarize, the cusp at fc = in the Airy 
PSF is not well fit by our default PSF characteriza- 
tion, leading to biases in the shape measurements. 
This can be remedied in a number of ways without 
invalidating our general approach. 

5.3. Accuracy for Asymmetric Galaxies 

The previous section discussed shape measure- 
ment accuracies with a well defined input shape 
and high degree of symmetry. However, real 
galaxy images do not have an unambiguous shape 
nor a perfect exponential radial profile. 

In this section, we measure the shapes of asym- 
metric galaxy images that have irregular isophotes 
(§4.2), for which there are no unambiguous defi- 
nition of shape. However, in weak lensing, it is 
not the accuracy of the individual shape measure- 
ment, but rather the ability to recover the shear 
from an ensemble of shapes, that is well-defined 
and needs to be accurate. In this sense, the ring 
test provides a measure of accuracy for assigning 
a shape for these irregular galaxies. 

In the ring test, the shape e is determined by 
GLFit measurements; the same image, before pix- 
elization, is rotated by a series of angles equally 
spaced between < /3 < 27r. Each rotated image 
is then sheared by 6+ = 0, 0.01 or 0.02, then pix- 
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Fig. 11. — Average shape (e+) obtained from a 
native fit ring test. The horizontal axis plots the 
measured galaxy ellipticity e. Each point is an av- 
erage from 10^ galaxies, while the error bar repre- 
sents the scatter in (e+); the dot-dashed line indi- 
cates the expected signal 1^5+. The (ex) response, 
not shown in this plot, exhibit a similar response 
as for the case 5+ — 0, since (5x = 0. 
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Fig. 12. — Calibration factor m in native fits (no 
PSF). The vertical axis plots the normalized shear 
signal {e+) /{TZS+) from Figm-e 11 (the data for 
S+ = 0.01 and 0.02 have been combined). The 
offset from unity is the calibration error. 



elized for the ring test. With 100 different orien- 
tations and 1, 000 (or more) realizations each, the 
uncertainty in the mean shape is as^ = cr^^ /VTO^, 
where (Te+ is the scatter in the e+ measurements 
due to photon statistics. The ring test has no 
intrinsic-shape noise. 

We test both the native and deconvolution fit 
cases at a series of different significance v, resolu- 
tion Tf), and galaxy ellipticity e. The asymmetric 
galaxies of different shape e are distinct combina- 
tions of exponential-disk and deVaucouleur pro- 
files; they are not sheared versions of the same 
underlying shape. 

Additionally, the shear calibration factor m and 
additive bias c can be measured with the ring test. 
The additive bias is the non-zero offset of the av- 
erage shape when setting 5 = in the ring test, 
and is typically present when epsp 0. The cali- 
bration factor m is then obtained by repeating the 
ring test with S ^ Q. 

5.3.1. Native Fits 

Figure 11 shows the mean shape (e+) as a func- 
tion of the galaxy ellipticity e for native fits. The 
mean shape is expected to follow the curve TZS-i-. 
The approximation TZ = 1 — e^/2, valid to 0{S^), 
is more than adequate here. For the various signif- 
icance {ly = 10, 20, or 50), the mean shape follows 
the expected curve very well. In particular, we 
see that there is no significant additive bias c seen 
in the 5+ = data points, to an uncertainty of 

Figure 12 plots the calibration factor m by nor- 
malizing (e+) for the 5+^0 cases to the expected 
signal 7?.5+, m = {e+) / {TZ5+). The calibration 
is well within 1% of unity aX v = 50, but at lower 
detection significance, there is a systematic under- 
estimation of the order m ~ 1 — Q.b/v. 

5.3.2. Deconvolution Fits 

The ring test with asymmetric galaxies is sim- 
ilarly performed with deconvolution fits, with a 
circular Airy PSF of epsF = and an elliptical 
PSF with epsF = 0.1. Figure 13 shows the addi- 
tive bias c in the shape average (e^) when 5+ ^ Q 
and epsF = 0.1, while Figure 14 show the calibra- 
tion factor m for (5+ = 0.02. Both Figures show 
various cases of (r^, v) pairs, which traces the 99% 
convergence contour from Figure 4. 
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Fig. 13. — Additive bias c due to PSF anisotropy. 
The average shape (e+) is obtained from a decon- 
volution fit ring test with an Airy PSF of epsF = 
0.1 for (5+ = 0, at various significance v and reso- 
lution Tf). The average shape is expected to he at 
zero (dot-dashed line); the offset from zero is the 
additive bias c. The dotted hues indicate ±1% of 
the PSF anisotropy. Each point is an average of 
lO'^ galaxies, (a) Additive bias at epgp = 0.10, 
for various (t^, r-h) pairs. (b) Additive bias at 
in = 1.2, 1/ = 40) for epsF = 0.05,0.10 and 0.15, 
with a scaling factor normahzed to epsp = 0.10. 
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Fig. 14. — Multiphcative bias in deconvolution 
fits, for various (y^r\^ pairs. The shear is ^+ = 
0.02, and the Airy PSF is circular, epsF,+ — 0. 
Each point is an average of 10^ galaxies. The 
dotted line above and below indicate ±1% of the 
shear. (Currently running more simulations to de- 
crease the error bar.) 

Figure 13 (a) shows c as a function of e at var- 
ious (rf), v) pairs. There is an e dependence to the 
additive bias, which increases in magnitude with 
decreasing v. Comparing with Figure 9, we see 
that the additive bias is essentially identical to the 
PSF bias found in symmetric galaxies, indicating 
that the PSF bias is independent of the galaxy 
symmetry, or the orientation of the galaxy shape 
over the ensemble average. Figure 13 (b) shows 
that for a given (rf,, v)^ the additive bias c scales 
linearly with epsF for reasonable epsp values. 

Figure 14 shows the calibration factor m for 
epsF = 0, 5+ = 0.02. At = 100 and better, the 
measured shear is within 1% of the input shear 
signal. As the significance decreases (y < 50), 
there is a multiplicative bias for galaxy shapes over 
e>0.6, whose magnitude is about —2.5% at e = 
0.8 (axis ratio a/b — 3). Similar results are found 
for m for the case epgp 7^ 0, where m = ((e+) — 
c)/7?.(5+, but with increased error bars due to the 
additional uncertainty in c. 
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6. Discussion 

6.1. Limitations of GLFit and Implication 

to WL Surveys 

There are limitations to GLFit with impHca- 
tions for WL surveys. If GLFit is to be used for 
shape analysis, the sampling rate must be such 
that the PSF minor-axis FWHM is sampled by 
> 2.8 pixels; otherwise, GLFit fails to converge to 
a unique shape (§5.1.1). 

Similarly, the PSF-smeared objects which pro- 
duce useful shapes must have a minimum reso- 
lution that depends on the significance (or vice- 
versa) in order to converge to a measurement of 
the pre-seeing galaxy shape. If we require shape 
measurements of all values of e to be successful, 
then Tf, > 1.8 is required at u = 20, or rt > 1.2 at 
1^ > 4:0 (§5.1.2). This implies that improvements 
in WL statistical accuracy are rapidly limited by 
the resolution of the galaxies: although going deep 
in the exposure increases the S/N of the detected 
objects and reduces the scatter and systematic in 
the shear estimate, the required S/N grows very 
rapidly for poorly resolved images, and it becomes 
difficult to produce additional useful WL data for 
poorly resolved galaxies. 

We have also found that larger, well resolved 
galaxies (r'b> 5) exhibit a problem with the decon- 
volution method when the PSF is an Airy func- 
tion, due to the poor approximation of the Airy- 
functions wings by the GL expansion. A simple 
remedy is available: smooth the PSF and galaxy 
images before measuring both. Another solution 
is to choose an alternative to Gauss-Laguerre as 
the PSF decomposition basis hmctions, one that 
is better suited to describe an Airy function. How- 
ever, a non-GL decomposition will make the de- 
convolution process excessively complex; another 
solution is to apodize the telescope to suppress the 
Airy wings. 

6.2. Comparison to Other Methods 

Kuijken (2006) offers an excellent description 
of the difference between various WL techniques; 
we refer the reader to this paper for the difference 
between, for example, the KSB method and Kui- 
jken's method, which is applicable to the difference 
bewteen KSB method and the EGL method. 

Our EGL method is similar to Kuijken's polar 



shapelet method (Kuijken 2006). What the meth- 
ods have in common are: the deconvolution of the 
PSF, which in principle allows for any PSF effects 
to be removed; forward fitting, which allows error 
propagation, and hence an error estimate to the 
measured shape; and the definition of shape as 
shear, which has a well-defined shear transforma- 
tion. All of these features contribute to a better 
shear accuracy. 

The differences between EGL and Kuijken's 
methods are subtle. The first difference is that 
Kuijken method works the deconvolution and 
shearing in shapelet-coefficient space. Our EGL 
method determines the shape-as-shear by itera- 
tively fitting within the pixel space. Secondly, 
Kuijken's method obtains the shape using a shear 
transformation valid to first-order in e, which can 
be off by up to 10% at e ~ 0.9 {g = 0.6). Our 
method uses basis functions that are elliptical, 
i.e. the shear transformation is valid to all or- 
ders, which allows for the shape to be measured 
accurately for any e. The third difference is that, 
in Kuijken's method, only the m = terms are 
used to describe the galaxy. This allows for less 
coefficients needed, and hence is efficient. In com- 
parison, our method obtains the full set of coef- 
ficients to the specified order, which theoretically 
is more accurate under shearing or deconvolution. 
Another difference is in the choice of scale ra- 
dius. Kuijken's method uses 1.3 times the best-fit 
Gaussian a for the basis function scale radius, 
which apparently is optimal for the first-order 
method; EGL uses the best-fit elliptical Gaussian 
size, which gives the optimal sensitivity to small 
shear (BJ02 §3.1). 

Unlike the KSB method, our EGL deconvolu- 
tion method, the Kuijken (2006) polar shapelet 
and the Kaiser (2000) methods do not rely on the 
approximation that the PSF is circular up to a 
small linear "smear." These methods are expected 
to do well on theoretical grounds; their differences 
(at least with Kuijken 2006) seem to matter at the 
1% level. 

6.3. Comparison to Other Tests 

As mentioned earlier, the main difference be- 
tween this paper and previous end-to-end tests 
(Erben et al. 2001; Bacon et al. 2001; STEP) is 
that we test the individual steps of the analy- 
sis. The difference is mainly in the distribution of 



20 



galaxies: these end-to-end tests mimic the distri- 
bution to an actual WL field image, while we test 
at specific S/N and galaxy size (resolution) sets. 
By controlling the exact distribution of galaxy 
shapes, we also eliminate shape noise in the shear 
estimate, and are able to quantify the effect that 
noise or resolution/sampling has on the shear. 

The main differences from other "dissection" 
methods (HS03; Kuijken 2006) are the inclu- 
sion of and quantifying noise effects, testing with 
asymmetric galaxies, and determining the lim- 
its of shape-measurement convergence. HS03 
demonstrates the calibration errors for a vari- 
ety of PSFs (Gaussian, seeing-limited, diffraction- 
limited), while we only test for Airy PSFs as the 
worst-case scenario. In HS03, the galaxy types 
were varied as well (Gaussian, Exponential, de- 
Vaucouleurs), but they are symmetric, which can 
mask potential problems in shape measurements. 
Their simulation focuses on the calibration accu- 
racy, but does not quantify the additive error or 
include noise effects. Kuijken (2006) includes an 
analysis with noise, but does not quantify them; a 
variety of Sersic indices (for galaxies) are tested, 
but are symmetric. Their PSFs do not test the 
diffraction-limited (Airy) case, although asymmet- 
ric PSFs, which are not considered in this paper, 
are tested. 

6.4. Accuracy with EGL Method 

With the EGL method, the systematic due to 
the shape measurement method is well under the 
0.1% level under ideal native- fit conditions, where 
the galaxy is symmetric and PSF is absent, as long 
as the minimum sampling and S /N criteria (minor 
axis width > 2.8 pixels, v > 10) are met. 

The deviation from this high accuracy comes 
from real-life complications. In native fits, when 
the galaxies are not elliptically symmetric, the 
shear calibration factor degrades by ^ —O.b/v as 
S/N is decreased. 

The presence of a PSF further complicates the 
matter; deconvolution with PSFs of truncated 
coefficients becomes a source of error for well- 
resolved objects with an Airy PSF (up to 0.5%), 
but solutions to this poor approximation of the 
PSF are straightforward. More generally, the 
circularly-symmetric Airy PSF itself introduces an 
e dependent systematic in the calibration factor. 



at worst ^ —3% at e = 0.8, ly = 20. A non- 
circular PSF adds a shear additive error that is 
proportional to the shape of the PSF, which also 
increases with decreasing S/N (up to ±3% of epsF 
a-tv = 20). For the deconvolution fits, the sign and 
magnitude of the calibration and additive error is 
approximately the same for both the symmetric 
and asymmetric galaxies. 

The systematics are typically inversely propor- 
tional to the S/N. At high S/N (t/ > 50), we have 
calibration error <1%, and the PSF anisotropy is 
> 99% suppressed. These errors become 4% 
problems as // 20. The systematics are a func- 
tion of the galaxy shape itself as well, with the 
systematic at low e having the reverse sign of that 
at high e. 

To attain 0.1% accuracy that is desired in fu- 
ture WL surveys, it is obvious that this shape mea- 
surement method itself needs to be refined. In ad- 
dition, the whole pipeline needs to be rid of any 
systematics that are unrelated to the shape mea- 
surement, such as crowding or selection effects. 

7. Conclusion 

The elliptical Gauss-Laguerre decomposition is 
one of the most stringently tested methods to char- 
acterize shapes of galaxies. With the EGL decom- 
position, shapes are measured without the need for 
empirical shear/smear polarizabilities, and PSFs 
arc removed by deconvolution. The shear, ob- 
tained from averaging an isotropic ensemble of 
galaxy shapes, is highly accurate due to the defi- 
nition of shape as shear. 

We have demonstrated that the EGL method 
allows shear recovery of unprecedented accuracy, 
and quantified its degradation due to PSF, trunca- 
tion of the EGL decomposition, image noise, and 
sampling rate/resolution. However, the current 
work is limited to the extraction of shapes; fur- 
ther work, including the full pipeline analysis, is 
required for attaining 0.1% accuracy in shear esti- 
mation. 
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